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Abstract.  The  immersed-boundary  method  can  be  used  to  simulate  flows 
around  complex  geometries  within  a Cartesian  grid.  This  method  has  been 
used  quite  extensively  in  low  Reynolds-number  flows,  and  is  now  being 
applied  to  turbulent  flows  more  frequently.  The  technique  will  be  discussed, 
and  three  applications  of  the  method  will  be  presented,  with  increasing 
complexity,  to  illustrate  the  potential  and  limitations  of  the  method,  and 
some  of  the  directions  for  future  work. 


1.  Introduction 

The  increase  in  computer  speed  achieved  over  the  last  few  years  has  made 
computational  fluid  dynamics  increasingly  useful  and  widespread  as  a tool 
to  analyze  and  design  flow  configuration.  Complex  geometries,  however, 
present  an  obstacle  even  to  present-day  computers,  since  the  use  of  body- 
fitted  meshes  (structured  or  unstructured)  significantly  increases  the  cost 
of  a calculation  in  terms  of  both  computational  speed  and  memory  require- 
ments. 

An  alternative  method  that  may  be  cost-efficient  in  many  situations  is 
the  “immersed-boundary”  method.  This  technique  is  based  on  the  intro- 
duction of  body  forces  distributed  throughout  the  flow  that  mimic  the  effect 
that  a solid  body  would  have  on  the  fluid.  This  approach  allows  the  use 
of  codes  in  Cartesian  coordinates,  which  present  significant  advantages,  in 
terms  of  speed,  accuracy  and  flexibility,  over  codes  that  employ  body  fitted 
grids. 

This  idea  has  been  widely  used  in  haemo-dynamics  and  bio-fluids  en- 
gineering: two-  and  three-dimensional  calculations  of  the  flow  in  the  heart 
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were  reported  by  Peskin  [13,  14]  and  McQueen  and  Peskin  [8,  9].  In  these 
calculations  the  motion  of  the  boundary  was  determined  by  the  fluid  it- 
self, so  that  the  boundary  had  to  be  modeled  as  a set  of  elements  linked 
by  springs.  In  cases  in  which  the  boundary  motion  is  known  a priori , the 
problem  can  be  significantly  simplified. 

Goldstein  et  al.  [4]  proposed  a feedback  forcing  mechanism  that  forces 
the  fluid  velocity  Ui  to  approach  the  velocity  of  the  solid  boundary,  Vj,  on 
the  boundary  itself.  Consider  the  incompressible  Navier-Stokes  equations: 


^1=  0 

Ox,  ’ 

dui  d 

dt  + dXi  {ujUi) 


1 dp 

pdxi 


+ vV2Ui  + fi . 


(1) 

(2) 


Goldstein  et  al.  [4]  assigned  a force  field 


= af  [ [ui(xSti,t)  -Vi(xSti,t)]dt'  + /3f[ui{xSti,t)  -Vi(xs 
Jo 


i>  ^)]  > 


(3) 

where  ay  and  fif  are  two  negative  constants,  and  xs^  are  the  coordinates 
of  the  solid  surface.  The  net  effect  of  this  force  is  to  tend  to  annihilate  the 
velocity  difference  ut  — Vi.  The  flow,  in  fact,  responds  to  the  forcing  as  a 
damped  oscillator  (see  [4]  and  [3]  for  an  in-depth  discussion  of  this  issue); 
the  frequency  of  the  oscillator  is  oc  |ay 11/2,  whereas  its  damping  coefficient 
is  oc  f3f/ \a\1!2.  This  implies  that,  in  order  to  enforce  the  no-slip  condition 
effectively,  ay  and  /3y  must  have  large  magnitudes  (larger  than  the  highest 
frequency  in  the  flow),  which  may  make  the  equations  stiff. 

Recently,  Mohd-Yusof  [11]  proposed  the  “direct  forcing  method,”  which 
assigns  a force  field  given  by 
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dx 


pdxi 


Vi  Ui 

At 


(4) 


(where  the  dependence  on  xSj  and  t has  been  omitted).  This  forcing  im- 
poses directly  the  desired  velocity  on  the  immersed  boundary,  and  has  the 
advantage  (over  the  feedback  forcing  method)  that  it  does  not  require  sig- 
nificant reductions  in  the  allowable  time-step.  It  was  extensively  tested  in  a 
staggered  finite-difference  code  by  Fadlun  et  al.  [3],  who  derived  an  interpo- 
lation scheme  to  be  used  when  the  boundary  does  not  coincide  with  a grid 
line.  Verzicco  et  al.  [15]  applied  this  method  to  the  large-eddy  simulation 
of  high  Reynolds  number  turbulent  flows  by  calculating  the  flow  inside  an 
IC  engine. 

In  the  present  paper  additional  applications  of  this  method  will  be  pre- 
sented, and  the  potential  and  limitations  of  the  technique,  as  well  as  issues 
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Computed  velocity  j_ 
| outside  the  body 

Interpolated  velocity 
....at  the  nearest  outside  point!. 


Figure  1.  Interpolation  method  used  to  apply  the  forcing. 


that  require  further  study,  will  be  discussed.  After  a brief  review  of  the 
governing  equations  and  of  the  numerical  scheme  used,  three  test  cases  will 
be  shown:  a low  Reynolds-number  flow  over  a cylinder  in  the  presence  of  a 
moving  surface  (Wannier  [16]  flow),  the  flow  over  a circular  cylinder  at  low 
Reynolds  number,  and  the  bypass  transition  on  a flat  plate  caused  by  the 
interaction  between  the  boundary  layer  and  the  wake  of  a circular  cylinder. 

2.  Problem  formulation 

Governing  the  flow  are  the  incompressible  Navier-Stokes  equations  (1-2). 
The  flow  solver  is  a standard  2nd-order  accurate  method  on  a staggered 
mesh  [1],  The  fractional  time-step  method  [2,  6]  is  used  and  a 2nd-order 
accurate  Adams-Bashforth  method  is  employed  for  the  time  advancement. 
A non- reflecting  boundary  condition  [12]  is  used  at  the  outflow,  and  periodic 
boundary  conditions  in  the  spanwise  direction.  The  inflow  and  free-stream 
conditions  depended  on  the  case  studied. 

The  direct  forcing  (4)  was  used.  Since  the  immersed  body  does  not 
follow  a grid  line  the  interpolation  method  proposed  by  Fadlun  et  al.  [3]  is 
used.  The  forcing  is  imposed  not  at  the  surface  itself,  but  at  the  first  point 
outside  it  (see  Fig.  1)  and  the  solid  body  velocity  in  (4)  is  replaced  with 
the  velocity  V)  obtained  by  a linear  interpolation  between  the  computed 
fluid  velocity  outside  the  body,  tq,  and  the  desired  body  velocity  V).  This 
method  has  several  advantages:  first,  it  has  been  shown  to  be  fully  second- 
order  accurate  in  time  [3];  therefore,  it  is  consistent  with  the  second-order 
differencing  scheme  used  by  the  solver;  secondly,  since  it  assumes  that  the 
velocity  profile  is  linear  near  the  body,  it  implies  homogeneous  Neumann 
boundary  conditions  for  the  pressure  (see  the  Appendix  in  [3]  for  a full 
discussion  of  this  issue).  This  last  feature  is  very  important  in  the  framework 
of  fractional  time-step  methods,  since  it  implies  that  the  corrector  step  does 
not  result  in  a modification  of  the  body  velocity  imposed  through  the  forcing 
in  the  Helmholtz  step.  On  the  other  hand,  the  assumption  that  the  velocity 
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Figure  2.  Wannier  flow  test  case  (Wannier  1950).  (a)  Computational  domain  and  com- 
puted streamlines;  (b)  L\  and  Li  norms  of  the  error,  e,  for  u and  v velocity  components. 
N is  the  total  number  of  grid  points. 


profile  is  linear  over  the  first  two  layers  of  cells  outside  the  immersed  body 
requires  the  use  of  a very  fine  mesh  in  the  vicinity  of  the  body. 

3.  Results  and  discussion 

3.1.  WANNIER  FLOW 

A straightforward  way  to  verify  the  accuracy  of  the  proposed  methodology 
is  to  compute  a flow  containing  a curved  immersed  boundary  for  which 
an  analytical  solution  exists.  The  case  considered  here  is  the  Stokes  flow 
around  a cylinder  in  the  vicinity  of  a moving  wall  (see  Fig.  2).  An  analytical 
solution  for  this  case  was  derived  by  Wannier  [16].  The  streamlines  for  this 
flow  are  shown  in  Fig.  2a.  Three  computations  on  gradually  finer  uniform 
grids  (32  x 32,  64  x 64,  and  128  x 128)  were  conducted.  The  L\  and  L2  norms 
of  the  error  (the  difference  between  the  computed  and  analytical  solution) 
are  shown  in  Fig.  2b  as  a function  of  the  total  number  of  points  N . The 
error  decreases  with  a —2  slope  indicating  that  the  proposed  methodology 
is  second  order  accurate. 

3.2.  FLOW  OVER  A CIRCULAR  CYLINDER 

The  next  test  case  examined  is  the  flow  over  a circular  cylinder  at  Re  d — 
UrXjD  jv  — 300  (where  D is  the  cylinder  diameter  and  the  free-stream 
velocity).  Two  calculations  will  be  compared:  a 2D  one  that  used  400x200 
points  in  the  xz— plane,  and  a 3D  one,  with  the  same  mesh  in  the  xz— plane, 
and  48  points  in  y.  The  computational  domain  was  60  x 27r  x 30,  and  the 
cylinder  center  was  at  xc  = 10,  zc  = 15  (all  lengths  are  made  dimensionless 
by  D,  all  velocities  by  Uoo).  The  grid  was  stretched  both  in  the  x—  and 
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Figure  3.  Flow  over  a circular  cylinder,  Rep  = 300.  Velocity  statistics,  (a)  U,  (b)  W, 
(c)  (uV),  (d)  ( u'w ').  Reference  data  from  Ref.  [10]. 


z-directions;  the  last  1/6  of  the  domain  (which  required  only  10  grid  points 
in  x ) formed  a sponge  region,  used  to  minimize  the  upstream  propagation  of 
disturbances  due  to  the  convective  outflow  conditions.  A uniform  velocity 
profile  was  imposed  at  the  inlet,  and  slip-wall  conditions  were  applied  at 
z = 0 and  z = 30. 


The  velocity  statistics  are  shown  in  Fig.  3.  Here  and  in  the  following  the 
angle  brackets  denote  averaging  in  time  and  in  the  spanwise  direction.  The 
3D  calculation  is  in  very  good  agreement  with  the  reference  data  by  Mittal 
and  Balachandar  [10].  At  this  Reynolds  number,  three-dimensionality  is 
observed  in  the  wake,  which  is  evidenced  in  the  visualization  in  Fig.  4, 
which  shows  iso-surfaces  of  the  second  invariant  of  the  velocity-gradient, 
tensor, 


Q = - 


1 din  duj 

2 dxj  dx{ 


2 Ojjfljj) , 


(5) 


(where  flij  is  the  anti-symmetric  part  of  the  velocity  gradient  tensor).  The 
condition  Q > 0 identifies  effectively  the  regions  of  coherent  vorticity  [5]. 

Figure  4 shows  the  formation  of  an  instability  on  the  initially  2D  rollers, 
and  the  presence  of  quasi-streamwise  rib  vortices  joining  the  rollers.  The 
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Figure  4 ■ Flow  over  a circular  cylinder,  ReD  = 300.  Isosurfaces  of  Q = 0.6. 


magnitude  of  the  spanwise  Reynolds  stresses  (v'v1).  in  this  calculation,  how- 
ever, remained  significantly  smaller  than  the  other  two  normal  components, 
which  allowed  the  2D  calculation  to  give  reasonable  results. 

The  effects  of  the  grid  resolution  near  the  obstacle  are  shown  in  Fig.  5. 
The  cell  Reynolds  number  is  defined  as 

Rec  = (A*2  + Af)  (v?  + m>) , (6) 

where  u and  w are  the  instantaneous  velocities,  and  Ax  and  Az  the  grid 
spacings.  If  the  mesh  is  insufficiently  fine  (Rec  > 30),  some  oscillations 
can  be  observed  that  initiate  along  lines  at  ±45°  on  the  cylinder  (they 
are  especially  visible  in  the  w contours,  Fig.  5b).  Refining  the  grid,  thus 
reducing  Rec  reduces  the  size  of  this  oscillation  (Figs.  7a  and  b).  Two- 
dimensional  interpolation  schemes  that  use  both  the  points  indicated  by 
the  diamonds  and  those  indicated  by  squares  in  Fig.  6 to  determine  Vi 
have  also  been  found  (Verzicco,  private  communication,  2001)  to  reduce 
the  amplitude  of  these  oscillations. 

3.3.  WAKE/BOUNDARY-LAYER  INTERACTION 

Wakes  interact  with  laminar  boundary  layers  in  many  applications  of  en- 
gineering interest,  for  example  on  the  leading  edge  of  multi-component 


-2-1  0 1 2 3 4 5 

.V 

Figure.  5.  Flow  over  a circular  cylinder,  Rep  = 300.  Coarse  (200  x 100)  2D  calculation, 
(a)  Contours  of  the  cell  Reynolds  number;  (b)  contours  of  w. 


Illfi 


Figure  6.  One-dimensional  vs.  two-dimensional  interpolation  stencils. 

airfoils,  or  inside  turbo-machincry.  The  interaction  of  the  turbulent  eddies 
present  in  the  wake  with  the  boundary  layer  itself  then  becomes  a primary 
driver  of  the  transition  process  in  the  boundary  layer  itself,  and  may  lead 
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Figure  8.  Sketch  of  the  wake/boundary-layer  configuration. 


to  transition  to  turbulence  at  fairly  low  Reynolds  numbers. 

The  configuration  examined  in  the  present  study  is  sketched  in  Fig.  8. 
A circular  cylinder,  with  its  axis  normal  to  the  stream  is  placed  above  a 
flat  plate.  The  cylinder  center  is  at  xc  = 10,  zc  = 3.2,  immediately  above 
the  leading-edge  of  the  plate,  which  was  also  at  xc  = 10.  As  in  the  previous 
case,  distances  are  normalized  by  D,  velocities  by  U0 c.  The  computational 
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domain  was  60  x 27r  x 20.  As  for  the  cylinder  calculation,  the  grid  was 
stretched  both  in  the  x—  and  z— directions  and  a sponge  region  was  used. 
The  Reynolds  number  based  on  cylinder  diameter  was  385.  The  configu- 
ration corresponds  to  Case  1 in  the  experimental  paper  by  Kyriakides  et 
al.  [7],  who  observed  significant  velocity  fluctuations  in  the  boundary  layer, 
starting  from  a point  approximately  six  diameters  downstream  of  the  cylin- 
der. These  fluctuations  are  generated  by  the  large-scale  convective  motion 
of  the  vortices,  and  do  not  die  down  after  the  wake  has  weakened,  but  de- 
velop into  a turbulent  boundary  layer  despite  the  fact  that  the  Reynolds 
number  is  very  low. 

The  distribution  of  the  streamwise  Reynolds  stress,  (u'u'),  as  a function 
of  x is  shown  in  Fig.  9.  A sudden  increase  of  (u'u1)  can  be  observed  to  begin 
at  x = 8,  indicating  the  beginning  of  transition.  This  result  is  consistent 
with  the  observations  of  Kyriakides  et  al.  [7],  who  defined  the  onset  of 
transition  as  “the  a:— location  where  the  velocity  signal  at  the  same  height 
above  the  plate  loses  its  sinusoidal  character”,  and  found  that  transition 
occurs  at  x = 7.4. 

The  velocity  profiles,  shown  in  Fig.  10a  at  several  locations,  initially 
resemble  a Blasius  profile  merging  into  a wake  near  the  cylinder.  As  the 
wake  widens  and  interacts  with  the  boundary  layer,  a logarithmic  layer 
begins  to  establish  itself,  indicative  of  transition  towards  turbulent  flow. 
This  transitional  behavior  is  also  observed  in  the  trace  of  the  Reynolds 
stress  tensor,  q 2 (equal  to  twice  the  turbulent  kinetic  energy),  which  in 
the  latter  sections  establishes  a turbulent-like  distribution,  with  a peak  of 
magnitude  7 — 8tw  at  z+  ~ 10  — 12.  It  should  be  noted  that  this  quasi- 
turbulent  state  is  achieved  at  very  low  Reynolds  number:  the  boundary- 
layer  thickness  <5  (defined  as  the  distance  above  the  plate  at  which  the  first 
maximum  of  the  velocity  profile  occurs)  is  approximately  50-70  wall  units. 

A visualization  of  the  flow  is  shown  in  Fig.  11.  The  structure  of  the 


Figure  9.  Wake/boundary-layer  interaction,  {u'u')  distribution  at  z — 0.18. 
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Figure  10.  Wake/boundary-layer  interaction.  Turbulent  statistics.  Top:  mean  velocity 
profile;  bottom:  q2  = (u'u'). 


cylinder  wake  is  similar  to  that  highlighted  in  Fig.  4,  with  strong  spanwise 
rollers  that  exhibit  3D  instabilities  and  eventually  break  up,  and  smaller 
quasi-streamwise  vortices  in  the  braid  region.  The  contours  of  streamwise 
velocity  fluctuation  u'  on  a plane  parallel  to  the  wall  show  significant  levels 
of  fluctuations,  especially  near  kinks  in  the  rollers  belonging  to  the  lower 
row.  Quasi-streamwise  streaks  are  formed  around  x = 15,  whose  spacing  is 
approximately  100  wall  units. 

4.  Conclusions 

The  immersed-boundary  technique  has  been  presented  and  discussed.  Il- 
lustrative results  from  three  simulations  indicate  the  potential  of  this  tech- 
nique, which  allows  the  calculation  of  flows  around  complex  geometries 
without  requiring  a body-fitted  grid. 

If  appropriate  interpolation  methods  are  used  when  the  body  does  not 
coincide  with  a grid  line,  the  method  is  second-order  accurate.  However, 
some  care  must  be  taken  in  the  discretization  of  the  flow  field,  especially  in 
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Figure  11.  Wake/boundary-layer  interaction.  Isosurfaces  of  Q = 0.4  and  contours  of  u 
in  the  z ~ 0.18  plane.  Top:  perspective  view;  bottom:  view  from  above. 

the  vicinity  of  the  body.  We  have  observed  the  development  of  numerical  os- 
cillations when  the  cell  Reynolds- number  exceeded  values  of  approximately 
30.  These  oscillations  did  not  grow  or  change  position  in  time,  and  their 
effect  on  the  flow  field  downstream  of  the  obstacle  was  limited  in  the  cases 
studied.  This  was  observed,  for  instance,  in  calculations  of  the  flow  around 
the  cylinder  at  Rep  = 3500  were  carried  out  in  which  the  grid  could  not  be 
sufficiently  refined  to  satisfy  the  cell-Reynolds-number  requirement.  How- 
ever, it  is  not  known  whether  these  oscillations  might  give  rise  to  instabili- 
ties in  other  geometries,  or  at  higher  Reynolds  numbers.  The  development 
of  more  accurate,  multi-dimensional  interpolation  schemes  might  be  bene- 
ficial in  this  respect.  The  use  of  multi-block  methods,  or  embedded  grids, 
could  also  alleviate  this  problem. 

If  these  numerical  schemes  can  be  overcome,  the  immersed  boundary 
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method  appears  to  be  a useful  tool  for  the  simulation  of  flows  in  complex 
geometries  at  moderate  or  high  Reynolds  numbers.  This  is  confirmed  by 
the  increasing  number  of  studies  using  this  method  that  are  appearing  in 
the  literature. 
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